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Abstract 


Many well known techniques exist for smoothing the periodogram estimate of power 
spectrum viz Bartlett’s and Windowing methods. ATS is a three step procedure 
for the same. First a small amount of local averaging (’A’) is carried out on the 
periodogram. Then a variance stabilizing transform is applied (’T’) and finally the 
result is smoothed (’S’) using a non-parametric regression procedure. This method 
works well even when data distributions are not necessarily Gaussian. Simulations 
were carried out on different data distributions using ATS and Bartlett’s method. 
Simulation results show that in a number of cases the proposed ATS method is 
superior to the Bartlett’s method. 



Contents 


1 Introduction 1 

1.1 Preliminaries 1 

1.2 Spectrum Estimation Methods 2 

1.3 Thesis Organisation 3 

2 Review 5 

2.1 Periodogram estimate of power spectrum 5 

2.1.1 Bias of periodogram estimate 6 

2.1.2 Variance of the periodogram estimate 6 

2.2 Smoothed Spectrum Estimators 7 

2.2.1 Bartlett’s method of averaging periodogram 8 

2.2.2 Windowing method 9 

3 The ATS method 11 

3.1 Basic Theory 11 

3.2 The Spectral Estimation problem 12 

3.2.1 Distribution of the estimate 12 

3.3 The averaging step ’A’ 14 

3.3.1 Bias and variance of the averaged estimates 14 

3.4 The variance stabilizing transform ’S’ 15 

3.5 The smoothing process ’S’ 17 

3.6 Choice of various parameters in the loess algorithm 19 

3.7 Scaled Exponential Variables 21 


1 



4 Simulation Results and Conclusions 23 

4.1 Generation of the Random Signal 23 

4.2 Simulations 24 

4.2.1 Simulations for the Stabilizing Transform 24 

4.2.2 Simulations for Periodogram Smoothing 25 

4.3 Conclusions 29 

5 Future Scope 42 


11 



List of Figures 


4.1 Power Spectrum of the AR2 Process y[n] 26 

4.2 Estimate Variance vs Frequency. (a=2) 27 

4.3 Estimate Variance vs Frequency.(a=4) 28 

4.4 Exponential White Noise Spectrum Estimate by ATS and Bartlett’s 

Methods. (a=2) 31 

4.5 Exponential White Noise Spectrum Estimate by ATS and Bartlett’s 

Method. (a=4) 32 

4.6 Exponential White Noise Spectrum Estimate by ATS and Bartlett’s 

Method. (a=8) 33 

,4.7 Spectrum Estimate of a Narrow Band Signal by ATS Method.(a=2) . 34 

4.8 Spectrum Estimate of a Narrow Band Signal by Bartlett’s Method. 

(Average of Two) 35 

4.9 Spectrum Estimate of a Narrow Band Signal by ATS Method. (a=4) . 36 

4.10 Spectrum Estimate of a Narrow Band Signal by Bartlett’s Method. 

(Average of Four) 37 

4.11 Spectrum Estimate of a Narrow Band Signal by ATS Method.(a=8) . 38 

4.12 Spectrum Estimate of A Narrow Band Signal by Bartlett’s Method. 

(Average of Eight) 39 

4.13 Spectrum Estimate of Sinusoids in Uniform White Noise by ATS 

Method. (a= 2) 40 

4.14 Spectrum Estimate of Sinusoids in Uniform White Noise by Bartlett’s 

Method. (Average of Two) 41 


ni 



Chapter 1 


Introduction 


1.1 P r eliminaries 

The need for power spectrum arises in a variety of contexts in communication sys- 
tems viz measurement of noise for design of optimal linear filter, measurement of 
narrowband signal in wideband noise, estimation of system parameters using an 
input noise excitation and measuring the output. Starting from the problem of es- 
tiniation of power spectrum we are often reduced to the problem of necessity for 
accurate measurement and fast computations. For exact measurement of power we 
would require an accurately measured infinite record length of the random signal 
which in turn necessitates infinite amounts of computation. These requirements 
being impractical an approximate measurement is necessary. Such an approximate 
measurement of power spectrum raises the issue of what record length of the random 
signal is required for a given accuracy, what computational approach should be used 
etc. Practical answers to these questions can be obtained by using certain results of 
statistical estimation theory. 

In our discussion we shall use two properties of spectral estimators viz bias and 
variance. Generally spealdng in case of parameter estimation if d is an estimate of 
a parameter a then 

Bias = a — e[Q:] 

Variance = e{a — e[d])^ 
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An unbiased estimator is one for which bias is zero. This means that the expected 
value of the estimator is equal to the true value of the parameter to be measured. 
A small value of variance suggests that the probability distribution function of the 
estimator is concentrated around its mean(true value in case of an unbiased estima- 
tor). If for an estimator the bias and variance both tend to zero as the record length 
increcises then such an estimator is called a consistent estimator [16]. 

Communication Systems are generally required to handle a variety of signals in 
the presence of noise. To a large .extent the design of such systems depends on 
the statistical properties of both the signal and noise. In most cases the noise can 
be represented (or approximated) as stationary Gaussian random process with zero 
mean. Certain applications require that signals and noises which are approximately 
stationary but not Gaussian be studied for their autocovariance or power spectrum. 
Although in a non-Gaussian framework autocovariance and spectrum are no longer 
the only relevent statistics, they are usually the most important ones [4]. However 
in either case measurement of power is of interest. 

We have only considered the case of measurement of spectra of individual signals 
or noises. The measurement of cross-spectrum is not discussed here. 


1.2 Spectrum Estimation Methods 

Broadly there are two methods for estimation of power spectrum. 

• Classical power spectrum estimation based on the Fourier analysis of the mea- 
sured random signal. 

• Parametric power spectrum estimation wherein a statistical model is dio- 
sen(AR, MA or ARMA) for the random signal and further analysis is carried 
out. 

The classical method applies the Fourier Transform on the random signal. The 
magnitude of the Fourier Transform squared divided by the record length is then 
a estimate of power. This is actually the same as the Fourier Transform of the 
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ACF of the random signal. Since we are considering only a finite record length of 
the random signal we are actually working with a windowed signal and a windowed 
ACF. Windowing means that we are assuming the signal and ACF to be zero outside 
the chosen window(The window size depends on the chosen record length). Often we 
have more knowledge about the random signal so that a more reasonable assumption 
than the signal or ACF as zero outside the chosen record length can be made [14]. 
This is done by choosing an appropriate model for the random signal. The problem 
then reduces to the estimation of the parameters of this model. This parametric 
spectrum estimation method is not used in our work. 

As shall be seen in chapter 2 the cleissical method does not give a consistent 
estimate. To make it a consistent estimate and to smooth the spectrum there exist 
well known techniques viz the Bartlett and Windowing methods. This work investi- 
gates the application of a new method (ATS method) to the same problem [8]. The 
estimate obtained from this method is compared with the one obtained with the 
Bartlett’s method. Since the method is supposed to work for non-Gaussian signals 
as well, simulations are carried out on non-Gaussian signals. 


1.3 Thesis Organisation 

The organisation of the remaining chapters is as given 

• As mentioned previously chapter 2 reviews the existing methods of periodogram 
smoothing viz Bartlett and Windowing method. These estimators are for their 
bias and variance. 

• Chapter 3 introduces the ATS method. This chapter describes in detail the 
method and the various steps involved. 

• Chapter 4 contains the simulations carried out using the ATS and Bartlett’s 
method. Comparisons between the two methods are made from the plots of 
the estimated power spectrum. The conclusions based on these plots are also 
listed in this chapter. 
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• Chapter 5 discusses the scope for further work using this method and improve- 
ments are suggested. 



Chapter 2 


Review 


In this chapter we will review the existing methods of smoothing the periodogram 
spectrum estimate. We will compare the advantages and disadvantages of these 
methods in terms of their bias and variance. 


2.1 Periodogram estimate of power spectrum 

For estimation of ACF and power spectrum of a random signal using time averages 
the signal should be ergodic. For an ergodic signal the time averages will approach 
the statistical average for increasing record length [17]. For any random signal 
there are two generally used estimators of ACF; one providing an unbiased estimate 
while the other a biased one. Both however provide a consistent asymtotically 
unbiased estimate. Unfortunately the Fourier Transform of this estimate is not 
a consistent estimate of the power spectrum. This is because the variance of the 
estimate does not tend to zero for increasing record length N of the random process. 
However smoothing this estimate does provide a consistent estimate of the power 
spectrum. For a discrete time real random process a:[n] let Cjjifm] be the estimate 
of the autocovariance and In{u) be the estimate of its power spectrum. Then for 
length 2M + 1 of the ACF Cxx[m] we have 

M 

ijv(w) = Cxj.[m]exp"-'‘^ (2.1) 

Tn=l 
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Corresponding to this if x[n] has a length N and Fourier Transform 

N 

X{exp^’^) = a;[«] exp"-''^ 

71=0 

it can be shown that [16] eqn 2.1 is the same as 

IsM = i I Jf (exp'"f I 

This is known as the periodogram estimate of power. 

2.1.1 Bias of periodogram estimate 

From eqns 2.1 and 2.3 we see that the periodogram estimate is equivalent to the 
Fourier Transform of the biased auto covariance estimator Ca;jc[m]. Taking the expec- 
tation of eq 2.1 we have 

Af 

= 13 e[cx®M]exp~-'‘^ 
m=l 

But from [16] 

£[cxx[m]] = ^ I m 1< Af 

Here is the autocovariance of the random process x[n]. Due to the finite 

limits of summation and the factor the expectation of is not equal to 

the Fourier Transform of Therefore the periodogram is a biased estimator 

of the actual power spectrum Pxx{<^)- Thus £[/iv(a?)] can be viewed as the Fourier 
Transform of the product of and <f>xx{‘m)- The equivalent frequency domain 

representation is 

£l/wH] = i jT (2.4) 

where Wj 5 (exp^") = is the Fourier Transform of 

2.1.2 Variance of the periodogram estimate 

We first consider the case when a:[n] is a white Guassian noise sequence. For this 
case it can be shown that [10] 

e(x[fc]a:[/]a;[m]a:[n]) = o^k = l and m = nor 


( 2 . 2 ) 

(2.3) 
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k — m and I = nor 
k = n and I = m 
= 0 otherwise 


Then the covariance of the estimate is given by [16] 


Cov[1nM,In{i^2)] = £^ x [( 


4 sm[(a;i +a?2)iV/2] 2 


iVsm[(a;i +a;2)/2] 


f + { 


sm[(cui - 0^2 )-^/ 2] .21 
jVsin[(a;i — a?2)/2] 


(2.5) 


Substituting wi = in eq 2.5 we get the variance of the periodogram in case of 


white gaussian noise as 

Var[/„M] = .11 + (^n (2.6) 

From eqn 2.6 we see that Var[/7\r(w)] is of the order of <r^. Also the variance 
does not approach zero as N increases. Thus the periodogram is not a consistent 
estimate of the power spectrum. From eqn 2.5 we see that the periodogram estimates 
seperated by integer multiples of ^ are uncorrelated. From this property and the 
non-consistency of the estimate we can conclude that as N increases the fluctuations 
in the periodogram become rapid. 

For the case where x{n] is a coloured noise we assume that it is derived from a 
white noise with unit variance by passing it through a linear filter. The magnitude 
squared response of this filter is Fa;ar(a;), where Pxx{^) is the true value of the power 
spectrum of coloured noise. Let xjv[n] be a AT point sequence of the coloured noise. 


If /jv(w) is the estimate of periodogram for a;jv[w] then [13] 


sm[(a?i - W 2 )N/ 2 ] , 2 
JVsm[(a>i — 0)2)72] 
(2.7) 


Here too we see that the frequencies which are integer multiples of ^ apart are 


uncorrelated. Also the estimate as in the case of white noise is not a consistent 


estimate. 


2.2 Smoothed Spectrum Estimators 

We shall now see two techniques to obtain a consistent estimate from the peri- 
odogram viz Bartlett and Windowing methods. 
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2.2.1 Bartlett’s method of averaging periodogram 

In this method a number of independent estimates are taken and averaged to obtain 
a reduction in the variance of the estimate. The origin of this method is attributed 
to Bartlett; hence the name. 

For a random process x[n] a sequence of length N is divided into K segments of 
length M each such that N = KM. Therefore 

= x[n + iM — M] 0 <= n <= M — 1; 1 <= i <= K 


The K periodogram estimates are computed as 

1 M—l 

= T7 I S ^^[njexp-^**^ 

n=0 

If the values of the autocovaxiance function (l>xx[^] is small for m > M then we can 
assume that the estimates of are independent of each other. Then the Bartlett’s 
estimate is defined as 

^ t=i 

Substituting eqn 2.4 of sec 2.1.1 for we have 




1 

27rM 



sm[(aj — 6)M/2] 
sm[(aj — 6)12] 


)^de 


Thus the expected value of the Bartlett estimate is equal to the convolution of 
the true spectrum with the Fourier Transform of the triangular window function 
corresponding to the M sample periodogram. Thus the Bartlett estimate is also a 
biased estimate of the spectrum. However since the estimate is obtained by averaging 
independent estimates, the variance is 


Var \ B , M ] = ^ VarlMu ,)] = 4 [ P „(..-)’][1 + ( 2 - 8 ) 


As the number of independent estimates increases the variance reduces and starts 
approaching zero. Thus the Bartlett estimate is a consistent estimate of power. 
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We have seen that the bias of the estimate is the convolution of the true spectrum 
with a spectral window. The width of the main lobe of the spectral window is 
inversely proportional to M. The bias increases with the increased width of the 
main lobe, consequently the spectral resolution reduces. For a fixed recod length as 
the number of periodograms increases the variance of the estimate reduces but due 
to a decrease in M the spectral resolution also reduces. Thus we see that there is a 
trade-off between the bias and the variance. 

2.2.2 Windowing method 

An alternative approach to smoothing the periodogram is by convolution with an 
appropriate spectral window. If Sxx{<jj) is the smoothed estimate and W is the 
appropriate spectral window then 

If w[m] = — f W(exp-^“) exp^‘^ du; 

2i7C tt 

is a finite sequence of length 2M — 1 then 

(M-l) 

Sxx{(*}) = Ca;j:(m)io(m) exp~'’‘*^ 

m=-(Af-l) 

Here Cxx[m] is the estimate of autocovariance. w[m] should be an even sequence 
so that Sxx{i^) is even and real when the input sequence x\n] is real. Also since 
the power spectrum is a non-negative quantity a sufficient although not necessary 
condition to ensure that 5a:®(a;) is non-negative is W(exp^") >= 0. This condition 
is satisfied by Bartlett’s window but not by Hamming or Hanning windows; so 
although the later provide a better frequency resolution they cannot be used [16]. 

Substituting for £[//^(^)] we get 
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If M is small compared to N then W(exp^‘^) will be approximately constant com- 
pared to Ws(exp^‘^). Therefore 

We see that the windowing method is also a biased estimate. Increasing the spec- 
tral window width increases the smoothing but reduces frequency resolution of the 
spectrum. The covariance and variance of this estimate are [13] 

If W is narrow compared to variations in Pix(w) then 

Var[S„(u.)] = (2.9) 

We see that as N increases the variance approaches zero. Thus the windowing 
method also gives a consistent estimate of power. 

Thus we see that these two smoothing techniques obtain a reduction in vari- 
ance of the periodogram estimate. Also though both these techniques are biased 
estimators, they provide an asymptotically unbiased and consistent estimate of the 
spectrum. In the next chapter we wil see in detail the new proposed method for 
periodogram smoothing. 



Chapter 3 


The ATS method 


3.1 Basic Theory 

ATS technique is an approach to fitting curves and surfaces to data using a non- 
parametric regression procedure. The probability distribution of the data to be re- 
gressed may not necessarily be Gaussian. Suppose Zi = (xi, y.) for Xi = (xn, x,' 2 , ...Xtv) 
f^ i=l to n are measurements of r-f 1 variables. The Xi are r-dimensional data and 
yi is some function of a:,-. Let y = g{x) be the functional relationship between Xi and 
y,-. For eg Xi may be a r-dimensional random sample and y,- its probability density 
function, y,- may be a binary response that depends on x,-.(In this case g{x) may be 
e(y/x)) For r=l and x,- spectral frequencies, y,- may be the periodogram estimate of 
a time series. (This is the application which we are going to investigate) [8]. 

In the ATS method ’A’ stands for averaging, ’T’ for transformation and ’S’ for 
smoothing. This method can be employed only for those data distributions for 
which a suitable variance stabilization transformation(’T’) exists. For the first step 
a suitable value of ’a’(number of points to be taken for averaging) is chosen. The 
average y,- is calculated by taking ’a’ values of yt whose Xfc are closest to x,-. The 
corresponding average x,- is also formed to get a set of data points of the form (x,-, y,). 
To identify the Xk closest to x,- a suitable metric is to be chosen over the space of x,-. 
The averaging can be done either by using overlapping or non-overlapping type of 
windows. Next a variance stabilizing transformation ’T’ is applied to the averaged 
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data to get another set of data (i,-,T(y,)) . In the final smoothing step ’S’ the 
{ti = T{yi)) are smoothed as a function of x,- using a non-parametric regression 
procedure. If required an inverse transformation can be applied to undo the effect 
of the ’T’ step and obtain an estimate of g{x). 


3.2 The Spectral Estimation problem 

To define the problem consider a real discrete time random process x]n\. We con- 
sider a windowed sequence of length N. Then the periodogram estimate of power 
spectrum as defined section 2.1 is 

= 4 I ^ xlnjexp-^’*^ P 

■''' n=0 

Here we estimate the periodogram at discrete values of frequencies for i=0 to 

[—]. This is in fact the Discrete Fourier Transform(DFT). Let y,- = /jv(w») be the 
estimates of the periodogram at frequencies x,- = a;,-. If y(a>) is the spectral density 
to be estimated then we shall try and estimate the same using the ATS procedure 
described before. The power spectrum can be estimated for any value of frequency 
using this procedure. Here though we are assuming that the power spectrum is a 
smooth function of frequency. 

3.2.1 Distribution of the estimate 

We shall show that for i = 1 to [y] the random variables are chi-squared 
distributed with degree of freedom 2 and for w,=0 and a;,=7r, have degree of 
freedom 1 [13]. In case of real x[n] we can rewrite the expression for periodogram 
estimate as 

n=0 n=0 

We first consider the case when x[n] is a white Gaussian noise sequence. 

e[x[n]] = 0 

e[x[n]^] = <7^ 
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^-1 

LetA{u)) = gfn]co5(a?n) 

n=0 

N-\ 

andB{ijj) = a;[n]sm(a;n) 

n=0 

Thene[A{u)]= £[B(a>)] =0 


For our problem since y,- are to be computed only at frequencies a?,- = ^ we consider 
only these frequencies. Therefore 


Var[A{o}i)] = e[A(w,)^] 


Also Cov\A{uJk)i A{u}i)] 
And Cov[B{u}}.)i 

Cov[A{u>k),B{(vi)] 


Onrin 

cos^-^) = 

a^N k = 0 or ^ 

0 k^l 
0 

0 for all k and I 



k = 



Since A{u}i) and J5(a;,) are linear functions of normally distributed random variables 
they are also normally distributed. Therefore 

A\u}i) _ 2A\u;i) B\ui) _ 2B^{oJi) 

Var[A{ui)] ~ Na^ yflr[B(a;,)] “ Na'^ 

are chi-squared distributed random variables with degree of freedom one (Xj). 
Therefore the sum jg chi-squared with degree of freedom two. For 

i = 0 or y B{ui) = 0 therefore t/,- = is chi-squared distributed with degree 

of freedom one. We now consider the case when a:[n] is not a white noise process. 
However a non- white process(coloured) can be obtained from a white one by suitable 
filtering operation. If z[n] is a white noise process, x[n] is the coloured noise process 
and /i[i] is the impulse response of the filter then 

OO 

a:[n] = — i] 

— OO 

If Pxx{t^) and Pzz{oj) are the power spectrums of a:[n] and z[n] respectively then 

=1 H(i^) p P,.(u) = 0^ I H(u) p 
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Consider a finite record of length N of the coloured noise process. 

= 4 I S x[n]exp“-'‘^” P 

n=0 

= 4 1 - «]exp"^‘^” p 

n k 

= ^ \ J 24 k-n] exp-'“(''-"> exp"^'"*' P 
=1 1 = 

• ^n{^) _ 

<7| I H(uj) (2 (r2 p^^(a;) cr| 

Since is X 2 for a;,- = i=l to [ j] and ATj^ for w.-O or|; therefore also has 
the same distribution. Thus we have seen that the periodogram estimates obtained 
are independent of each other and are distributed. The frequendes other than 
w = 0 and w = | having degrees of freedom 2 are exponentially distributed random 
variables. For the ATS operation we drop the end frequendes 0 and | and consider 
only the other frequendes. 


3.3 The averaging step ’A’ 

In the averaging step we average the periodogram in ’b’ blocks of ’a’ consecutive 
Fourier frequencies. We use averaging with disjoint windows and since all quanties 
are unidimensional the metric used is arithmetic difference of the quantities under 
consideration. The average yi for the block is associated with the averaged 
frequency x/. The number of blocks thus obtained is 6 = Since all the 

quantities averaged are exponential random variables their average will be gamma 
distributed. Averaging with disjoint windows maintains the independence of the 
estimates. 

3.3.1 Bias and variance of the averaged estimates 

The averaging done in the ’A’ step is equivalent to the windowing technique of peri- 
odogram smoothing. Here a rectangular window(in the frequency domain) is used. 
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This leads to reduction in variance of the estimate, the reduction being proprtional 
to the value of ’a’. Thus the expressions of bias and variance derived for the window- 
ing step are valid for this. However since we are computing the periodogram only at 
discrete values of frequency the integration operation is replaced by the summation 
operation. Also the bias and variance are proportional to the frequency a>,-. The 
expressions for mean and variance as functions of the frequency oji are [15] 


mean{u)i) 


27r 




2 


E 




iy(exp^®) 


(3.1) 


VariiJi) = 


Plx{‘^i) r 1 

N ^TT 


(g-i) 

2 

E 






(3.2) 


Here W(exp^®l) = ^ for all 6. Prom the above expressions we see that the mean and 
variance depend on the frequency a;,-. Also the variance is related to the mean by 
the relation V’ar(a?i) = Kmeari^iijJi) for a suitable constant K. Since the mean and 
variance depend on the power content at that frequency, for a non-white process the 
variance will change with the frequency. Regression analysis is used to fit a curve 
of the form y,- = g{xi) -f where e,- are i.i.d random variables. Thus it is necessary 
to investigate whether there exists a transformation of some kind which will make 
the variance of the estimate a constant at all frequencies which is independent of 
the frequency and the power content at that frequency. 


3.4 The variance stabilizing transform ’S’ 

The purpose of this transformation is to change the scale of the observed data to 
make the data analysis more valid. The transformed variate should also have a con- 
stant variance. Ideally the transformed variable must have the following properties 

• Variance of the transformed variable must be unaffected by changes in the 
mean level. 


• The transformed variable should be normally distributed. 
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• The transformed scale should be one in which the arithmetic average is an 
efficient estimate of the mean. 


• The real effects should be linear and additive. 


The required transformation function depends on the mean, variance and distri- 
bution of the random variable. Consider a random variable x whose mean £{x) is 
a real variable having the range R. If is the variance of x, let it be related to 
the mean m through the relation cr* = /(m). Thus any fluctuations in the mean 
are going to affect the variance. In case of averaged periodogram the mean being a 
function of the power at that frequency, the variance is also then a function of power. 
In order to stabilize the variance ie make it a constant irrespective of the power or 
frequency it is necessary to transform the variable x through a transformation <(x) 
such that the variance of the transformed variable is independent of the mean. This 
means that we should achieve 


dy dcrl 

= 0 and = 0 
am dm 

for all m in R. Also starting from ^ = t'{x) it can be deduced by a summation 
procedure that Cy = t'{m)cr{m). This derivation is not mathematically sound and 
can be used if on application it gives satisfactory results. For our case we shall use 
this result to derive the necessary transformation function. 

After application of ’A’ step the relation between mean and variance is 


cr(m) = Km 


where K is a constant. If j/ = t{x) is the transformation being used then it is required 
that Var[y) = <ry=constant(say c). Thus we have 




ie 


c = - — Km 
dm 

dt = K'^^ ort — K'logim) 
m 
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In general we have t{x) = Klog{x). Thus for a chi-squared periodogram estimate 
which is smoothed by a window function the variance stabilizing transform is the 
logarithm function. More details on the logarithm transformation can be found in 
[9] [3] . As stated before the transformation function depends on the distribution of 
the random variable. In sec 3.1 we had stated that one of the other applications 
of this method is the case when y is a binary response on x. Here the required 
transformation is the inverse sine function [12]. 

3.5 The smoothing process ’S’ 

The final step consists of smoothing by regression. We have the transformed vari- 
ables U = log{yi) which are i.i.d random variables. We use a regression tedmique 
known as locally weighted regression(loess) [6]. 

Locally weighted regression or loess is a method of estimating a regression surface 
through a multivariate smoothing procedure. It consists of fitting a function of 
independent variables locally and in a moving fashion analogous to a moving average 
filter. Loess is different from the usual polynomial regression. Loess can be used 
to estimate a wider class of regression surfaces than the polynomial method. Since 
the spectum smoothing is an univariate problem we shall discuss only the univariate 
loess method. The same technique has been generalized and used for the multivariate 
case [7]. 

Consider a set of points of the form By applying locally weighted re- 
gression we shall estimate a new set of points This new set of points y,- 

are the fitted values at x,-. If the regression surface to be estimated is of the form 
Vi = gi^i) + £,■ then yi = g{xi). 

In order to describe the smoothing procedure we consider a window function 
having the following properties 

1. W(x) > 0 I X |< 1 

2. W{-x) = W(x) 

3. W (x) is a non-increasing function for x >= 0 
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4. Wix) = 0\x\>^1. 

Let / be a parameter(0 < / <= 1) used to control the amount of smoothing in the 
regression process. If U is the number of points of the form which are to be 

smoothed; let r = t// be rounded off to the nearest integer. For each Xj the weights 
Wk{xi) are defined for k=l to U. This is done by centering the weight function W at 
Xj and scaling it so that the point at which it first becomes zero is the nearest 
neighbour of x,-. Then j/,- is determined as the value of a degree polynomial fit to 
the data using weighted least squares with weights Wk(xi). This procedure is carried 
out for every x,- to obtain the values yf. After this a set of residuals Si = yi — yi is 
defined for each x^. New fitted values are computed using the same procedure just 
described but using a new set of weights 5»iojfe(x,). This procedure is carried out for 
a specified number of iterations after whicL a smooth fit is obtained. 

The assumption of smoothness allows us to use the points in the neighbourhood 
of (x,-,y,) to determine y,-. The values of the weight function Wk(xi) decreases as 
the distance from Xi increases. Thus the points which are close to Xj play a larger 
role in determining the value of y^ than those points which are farther away from 
X,-. Increasing the value of ’f’ increases the neighbourhood of influential points and 
therefore increases the smoothness of the fit. 

We shall now breifly describe the algorithm used in the loess procedure. For 
each X,-, h,- is the distance of the nearest neighbour of x,- ie hi is the smallest 
value of 1 Xi — x*. | for k=l to $. The weighting function is defined as 

Wk{Xi) = 

hi 

The following sequence of operations is carried out 

1. At the point (x,-,yi) we are fitting a order polynomial to obtain the 
value of y,- as y,- = The estimates of the polynomial coeffi- 

A 

cients ^j{xi) are determined using the method of weighted least squares. 
The ^s' are obtained by minimizing the value of 
S\i 

5D ^ki^diVk -fio- 
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Then 

y< = S 

3=0 

This is carried out for all Xj. 

2. Let B be a bisquare weight function which is defined by 

B(x) = (1 - x2)2 1 1 |< 1 
= 0 I X |>= 1 

If Si = yi — yi are the residuals from the fitted values; let ’p’be the median 
of these fitted values 1 e, |. We define robust weights as Sk = B{^). 

3. The steps (1) and (2) are repeated using the weights ^fciOfc(xj). Thus a 
new set of smoothed points y^ is obtained. This procedure is carried out 
for a certain number of iterations. 

The weight function is derived from the tricube function 

W(x) = (1- I X I X |< 1 
= 0 I X |>= 1 

The above procedure can be used to determine the value of y at any value of x(ie 
other than x,). Thus a smooth and continuous curve of the power spectrum can be 
obtained. 


3.6 Choice of various parameters in the loess al- 
gorithm 

There are four parameters which are to be considered in the algorithm. They are 


1. The order of the polynomial fit/d'. 

2. The weighting function to be used.'t/;'. 

3. The number of iterations, 'i'. 
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4. The smoothing parameter//'. 

The choice of 'd' is a trade-off between the need to reproduce variations in the 
spectrum and computational ease. The simulations done by us are performed for 
d = 2. However the algorithm can be used for any value of ’d’. 

The requirements of the weighting function 'w' have been stated before. The 
reasons are stated now. The first requirement is chosen because negative weights do 
not make any sense. The second requirement W{—x) = kr(x) is chosen to give the 
same importance to points which are at the same distance on both sides of Xj. W^(x) 
should be non-increasing with increasing x so that a point nearer to Xj will have 
more weight than a point which is farther from x,-. In addition Vr(x) should decrease 
smoothly from x = 0 to x = 1. The tricube function is chosen because it enhances 
a chi-squared distributional approximation of the estimate of error variance [6]. In 
our simulations we have used two iterations. It has been studied that two iterations 
give a good, robust fit. However some type of convergence criterion can be defined 
for the loess algorithm. Execution of the algorithm is carried out till this criteion is 
met. This method however has not been formulated or implemented. 

It had been stated earlier that this method can be used only for cases where the 
spectrum is a smooth function frequency. However the extent of smoothness of the 
curve is not known and no procedure has been adopted to measure it. Broadly it can 
be stated that the amount of smoothing must be enough to remove the fluctuations 
in the spectrum due to the statistical nature of the data and finite record length. 
But it should not be excessive so as to distort or smoothen out the actual variations 
in the spectrum. For eg for a narrow band signal whose power is concentrated in 
a narrow band of frequencies smoothing may be dangerous for the spectrum. The 
simulations were therefore carried out for different smoothing factors and the plots 
compared. 
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3.7 Scaled Exponential Variables 


To illustrate the problem of smoothing a periodogram consider the problem of 
smoothing scaled exponential variables [8]. Let 


Vi = PiWi 


-(n-l) 

2 2 


where Wi are independent standard exponentials and m are smooth functions of i. 
Let Hi = exp(Q; + i/3) 

We assume that the objective is to estimate the centre value of yu ie exp(a). If 
we assume that the record length n is of the form n = ab then we can divide the 
data in ’b’ disjoint blocks and apply the ’A’ step that averages the ’a’ values within 
each block. In the ’T’ step we take the logarithms and in the ’S’ step we shall 
perform smoothing by averaging by averaging these ’b’ quantities. This process 
should give a good estimate of a. Since exp(a + i^) is symmetric about the centre 
value, smoothing by averaging gives the same estimate as obtained by a linear fit. 
After the ’A’ and ’T’ steps we get 

(e-l) 

j/j = a + Paj + ln{f) + ln{ 2Wjk/a) + Zj 

where Wjk = Wj+k,fk = exp*^,/ = E ^.-u exp*^,</jfe = ^,Vjk = Zj = 

^”(1 + J2k9kVjk)- 

The value of j varies from 1 to ’b’. The ’S’ step gives y = ave{yj). We shall now 
derive the mean and variance of y and yj. Since Wi are independent exponentials 
E,._ E(a-i) W is gamma(a) distributed. Then [1] 

K—~ o ' 


e[ln{Y^W)] = xj;{a) 


Var[ln{J2W)] = ip'{a) 

k 

is the derivative of 0(a) and gamma’{a) is the derivative 
of gamma(a). To calculate the moments of Zj we expand it by its Taylor’s series 
expansion. 






r=l 
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Using the independence of V and ^ W we have 


^ <Ewy 


Consider the expression 


=n ^ • L r-0 ”• 


) 


= £(exp(^5;^l^)) = ni(l - eg)-^ = 

k r 

Here is the homogeneous symmetric polynomial in the g's which can be 
expressed in terms of the power sums Sr = I2r [H]- Also we have from the above 
expressions 

e(£gWY = r\hr 

k 

From the tables in [11] we get hi = Si = 0,2^2 = ■82,6^3 = 253,24^4 = 6S4 + Ssj. 
Thus we get 

-/ 7 'v _ , 2^3 _ 3(2s 4 + s^) 

roi "r o rii fyii "r ••• 


2a[2] ^ 3aW 


oW 


2\ ^2 2s3 11(2S4 + 52) 


4aW 


+ ... 


= a{a + l)(a + 2). ..(a + r — 1) 
From these we get the moments of y, and y = ave{yj) as 

£(y) = a + ln{f) + ^(o) - ln{a) + e(Z) 


nVar{y) = a‘tl?'{a) + aVar{Z) 

It turns out that for small values of ’a’ (a < 8 approx) the contribution of 6 [Z) and 
Var(Z) to the mean and variance of y is negligible. Thus the mean and variance 
of y mainly depend on ^-nd ’a’ itself (for small values of ’a’). Also the 

variance of tj ie the yj is approximately 

Var{tj) = 0^(a) 

Thus we see that the variance of the transformed variables tj is a constant which 
is independent of their distributions. These stabilized variances have been verified 
through simulations. 



Chapter 4 


Simulation Results and 
Conclusions 


This chapter contains the simulations carried out and the conclusions inferred from 
these simulations. 


4.1 Generation of the Random Signal 

The first task is the generation of the necessary random signal (data). For this 
the pseudo random number generator algorithm of Park and Miller was used [18]. 
This algorithm generates a pseudo random number which is uniformly distributed 
between 1 and 2®^ — 2. This number is then divided by 2^^ — 1 to generate a random 
number which is uniformly distributed between 0 and 1. (exclusive of the end-points) 
This uniformly distributed number can then be used to generate a wide range of 
distributions. For this purpose we make use of the following theorem 

If U is uniformly distributed in [0,1] then F~{U) is a random variable whose 
cumulative distribution function is F. 

Using this we can generate different distributions provided their cumulative dis- 
tribution function is known and is invertible. Thus in order to generate an exponen- 
tial random variable whose probability distribution function is f{x) = Aexp(— Ax) 
the necessary function F~{U) is ^ln{l — U). However since U is uniformly dis- 
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tributed so is (1 — {/). Therefore the required exponential random variable can be 
generated using The data to be generated has to be initialised with a 

seed value which is a negative integer. The average of the entire record gives an 
estimate of the mean which can be substracted from each value to get a zero mean 
random variable. For a zero mean random variable the autocovariance is the same as 
autocorrelation. For different seeds the data stream starts from different numbers. 
Similarly other types of data can also be generated. 


4.2 Simulations 

4.2.1 Simulations for the Stabilizing Transform 

In sec 3.4 we had learnt about the variance stabilizing transformation. We had 
seen that the ATS method can only be used if such a transformation exists. Also 
we had seen that for the periodogram the required transformation is the log func- 
tion. We shall verify through simulations if this transformation results in a constant 
variance that it claims. For this purpose it is necessary to perform the simulations 
independently a large number of times and the variance must be calculated from 
these simulations. Here we are assuming that the time averages will give us a good 
estimate of the variance. The following experiment was carried out 

1. For i=l to M do upto step 6 

2. Using seed value as -i generate uniform random numbers x[n] in [0,1] for 
n=l to N. 

3. Generate y[n] = — ^] + 2 ;[n] for n = 1 to N. 

4. Calculate the FFT of y[n] at frequencies w* = ^ for = 0 to y and 

estimate the power spectrum from this using eq.2.3. 

5. Divide the entire frequency range into ’b’ blocks of ’a’ frequencies each. 

I _ tjjg average value of the estimate and the average 

value of the frequency. Let /jv denote the estimate, then store this result 
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in which indicates that this is the power spectrum estimate for 

the frequency for the seed value, j varies from 1 to b. 

6. Take the'logarithm of /ivL/j*] and store it in TatL?,*]- 

Then the variance before the transformation at frequency j is calculated as 

1 1 A/ 

Similarly the variance after the log transform <tjj can also be calculated by using 
in place of /a/[7i*] the above formula. The above simulations are carried 
out for diflferent values of averaging ’a’. In our simulations the record lenth N was 
taken to be 1024. The variance calculations were done by averaging over 10000 
independent iterations. The data y[n] was generated from a uniformly distributed 
data a:[n] using a second order AR process with oi = 0.0 and 02 = 0.9025. This 
data has power spectrum as shown in fig.4.1. 

Fig.4.2 is the plot of variance versus frequency for averaging value a=2. Fig.4.2(a) 
shows the variance before the transform while fig.4.2(b) shows the figure after the 
transform. We can see that variance before the transform has the same profile as 
that of the power spectrum(see fig.4.1) while the variance after the transform is 
approximately constant over the entire frequency range. This constant variance is 
equal to 0.64 which agrees with the theoritical value 

^'(a) = 0.6445 at a = 2 

Fig.4.3 shows the plots for a=4. Here the variance after the transform stabilizes at 
0.28. Theoritically we have 

V»'(a) = 0.28375 at a = A 

4.2.2 Simulations for Periodogram Smoothing 

The next simulations are for the smoothing of the periodogram estimate. The 
smoothing was done using the ATS and the Bartlett s method. Here too the op- 
erations carried out are the same till the log transform step. After that the loess 
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Figure 4.1: Power Spectrum of the AR2 Process y[n] 

algorithm is used for smoothing the spectrum. The variables in the loess algorithm 
are the smoothing factor ’f ’ , degree of the polynomial to be fit ’d’ and the number 
of iterations ’t’. In order to estimate the spectrum by Bartlett’s method we divide 
the record length N in ’a’ disjoint records so as to have ^ as the length of each. The 
pefiodogram estimate of each block is estimated independently and the average of 
the ’a’ periodograms is the Bartlett’s estimate. Here ’a’ is the averaging factor used 
in ATS. We have used a record length N of 4096 data points. The degree of the 
polynomial fitted is 2 and the number of iterations carried out also is 2. 

Figs.4.4 to 4.6 are the plots for estimating the spectrum of exponential white 
noise.(with mean arrival rate A = 1) Fig.4.4 is for ’a’=2. In figs.4.4 (a), (b) and 
(c) the smoothing factors used are 0.01, 0.05 and 0.1 respectively. Fig.4.4 (d) is the 
Bartlett’s estimate as the average of two disjoint periodograms. Figs.4.5 and 4.6 are 
the plots of the simulations carried out for ’a’=4 and ’a’=8 respectively. Next the 
simulations were carried out on a narrow-band signal whose power spectrum is as 
shown in fig.4.1. The data was generated using the iterative procedure 

y[n] = — 0.9025y[n — 2] -f x[n] 

where a:[n] is uniform in [0,1]. Figs.4.7 and 4.8 are the plots for a=2. Fig.4.7 are the 
plots of the estimate using the ATS method for four smoothing factors of 0.01, 0.05, 
n 1 anrl fl ^ frmrth Qmnnt.hirnr fartor has been included to show how excess 
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Figure 4.3: Estimate Variance vs Frequency. (a=4) 
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smoothing can distort the spectrum. Fig.4.8(a) is the plot of the spectrum after 
the A and T steps so as to compare it with the estimates obtained after different 
values of smoothing. Fig.4.8{b) is the Bartlett’s estimate. Figs.4.9 and 4.10 are 
the plots for a=4 while figs.4.11 and 4.12 are for a=8. The first case of estimating 
the spectrum of white noise was a smooth spectrum while the second spectrum had 
some variations in it. In the third case we shall take a spectrum having sharp peaks 
or discontinuities in it. This gives us a fairly good idea of how the method works 
for a smooth spectrum at one end to a spectrum having sharp peaks at the other 
end. Therefore the third case we estimate the spectrum of sinusoids buried in white 
noise. 

x[n] = sm(a;in) + sin{ci} 2 n) + u[n] 
where t;[n] is uniform in [-0.5.0.5]. 

a>i = 0.5 and wj = 0.6 

In the plots the x-axis has the normalized frequency. Therefore Wi = 0.5 and 
0)2 = 0.6 correspond 0.0796 and 0.0955 on the X-axis respectively. Fig.4.13 shows the 
plots using the ATS method using smoothing factors 0.01, 0.05 and 0.1 respectively. 
Here we see that the smoothing is not able to resolve the two sinusoids. In fact 
fig.4. 14(a) which has the spectrum after the ’A’ and ’T’ steps the two sinusoids 
were clearly resolved. Smoothing thus has distorted the spectrum. The Bartlett’s 
method fig.4. 14(b) gives a good estimate. The method can also be tried on Gaussian 
distributions. 


4.3 Conclusions 

We had seen before that the ATS method can only be used in cases where a variance 
stabilization transformation exists. The importance of this transform has been veri- 
fied through simulations. Due to the transform the estimates have the same variance 
at all frequencies and hence a regression procedure can be used. Regression is used 
to fit a model of the type yi = si^*) "i” i.i.d. and have a constant 

variance. The log transform effects this constant variance. 
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As seen from the various plots the ATS method is quite efficient when the power 
spectrum is a smooth function of frequency. As seen in the case of the exponential 
white noise this method gives a better estimate than the Bartlett’s method. The 
Bartlett’s method estimate has a lot of fluctuations due to the statistical nature 
of the data and due to the finite record length. Regression helps to smooth these 
fluctuations. In the case of the narrowband spectrum the ATS method gives better 
results for small values of the smoothing factor ’f’. At these values the estimate is 
better than the Bartlett’s estimate. However for large values of smoothing factor it 
distorts the spectrum. Compared to this the Bartlett method gives better results 
in the sense that it does not distort the actual variations in the spectrum. However 
if the spectrum has sharp peaks then Bartlett’s method can be preferred over the 
ATS method. 

Thus we can conclude that the ATS method gives good results(better than 
Bartlett’s method) if the spectrum to be estimated is a smooth function of fre- 


quency. 
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Figure 4.5: Exponential White Noise Spectrum Estimate by ATS and Bartlett s 
Method. (a=4) 
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Figure 4.6: Exponential White Noise Spectrum Estimate by ATS and Bartlett’s 
Method. (a=8) 




Figure 4.7; Spectrum 


Estimate of a Narrow Band Signal by ATS Method. (a-2) 
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Figure 4.8: Spectrum Estimate of a Narrow Band Signal by Bartlett’s Method 
{Average of Two) 
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Figure 4.10: Spectrum Estimate of a Narrow Band Signal by Bartlett’s Method. 
(Average of Four) 
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Figure 4.11: Spectrum Estimate of a Narrow Band Signal by ATS Method.(a— 8) 
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Figure 4.12: Spectrum Estimate of A Narrow Band Signal by Bartlett’s Method. 
(Average of Eight) 
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Figure 4.14: Spectrum Estimate of Sinusoids in Uniform White Noise by Bartlett’s 
Method. (Average of Two) 



Chapter 5 


Future Scope 


The ATS method essentially consists of estimating a regression surface(or curve) by 
smoothing. Much research has been carried out and there exists lot of literature in 
this field. 

One important property of the ATS method (which has not been made use of 
in this work) is that the residuals obtained by regression can be used to determine 
thye adequacy of the fit. For this purpose the residuals are to be compared with the 
fitted values in the smoothed spectrum. Larger the smoothing factor larger will be 
the size of the residuals. Thus the residuals can be used as an important factor for 
deciding upon the amount of smoothing. 

Another method to choose the optimum value of ’f; is to find it iteratively. The 
method is as described. Let y«(/) be the fitted value at Xi for a given value of f 
with yi not included in the computations. The initial fitted value of / = /o is chosen 
by minimizing 

'^{yk-yk{f)Y 

k=l 

where n is the number of points in the regression. Once this value /o is found out 
the robustness weights Sk can be obtained as done in the loess algorithm. The next 
value of f is then chosen by minimizing 

t^k{yk-yMf)? 

This procedure is repeated till convergence takes place at the final value of f. This 
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technique has been used in some applications [6] although not for the periodogram. 

The loess procedure does local fitting using the method of least squares. Another 
method is spline smoothing for the estimation of regression surface [19]. Here a 
smoothing procedure using splines for the log periodogram is discussed. Essentially 
we have looked into two aspects. One is the choice of the optimal smoothing factor 
’f’. The other is an alternative to the locally weighted regression used in the ATS 
method. 
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